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In multicellular organisms, several cell states coexist. For determining each cell type, 
cell-cell interactions are often essential, in addition to intracellular gene expression dynam- 
ics. Based on dynamical systems theory, we propose a mechanism for cell differentiation 
with regulation of populations of each cell type by taking simple cell models with gene 
expression dynamics. By incorporating several interaction kinetics, we found that the cell 
models with a single intracellular positive-feedback loop exhibit a cell fate switching, with 
a change in the total number of cells. The number of a given cell type or the population 
ratio of each cell type is preserved against the change in the total number of cells, de- 
pending on the form of cell-cell interaction. The differentiation is a result of bifurcation of 
cell states via the intercellular interactions, while the population regulation is explained 
by self-consistent determination of the bifurcation parameter through cell-cell interactions. 
The relevance of this mechanism to development and differentiation in several multicellular 
systems is discussed. 
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1 Introduction 



Complex gene regulatory or protein networks are responsible for determining cellu- 



speci fic network structures called network motifs (|Shen-Orr et al 



discussed in the light of 


. 2002: 


Milo et al.. 


2002, 



2004 ) . Besides such motifs, several simple network modules are also considered to op- 



erate to g ive specific dynamical properties such as b i stabi 



behavior ( Ferrell Jr. and Machleder 



1998; 



Sha et al 



2003 



i ty, adaptation 



Tyson et al 



, or o scillatory 



2003J). Recent 



experimental results also suggest that such modul es provide a basis for cell differen t iation 



as studied in competence state in Bacillus subtilis (jSiiel et al 



2006 



Maamar et al 



20071 ) 



In multicellular organisms, several cell states coexist. Morphogenesis with differen- 
tiation into distinct cell types, however, is not an event of independent single-cellular 
dynamics, but occurs as a result of an ensemble of interacting cells. For determining 
each cell type, cell-cell interactions are often essential besides intra-cellular dynamics 
by functional modules at a single cell level. In fact, gene regulatory networks respon- 
sible for the early developmental process or the cell s pecification process of several kinds 



of organisms include many intercellular interactio ns 



2007 



E. H. Davidson et al. 



2002: 



Imai et al. 



2006 



Ben-Tabou de-Leon and Davidson 



Loose and Patient 



2004 



Swiers et al 



20061 ). The importance of cell-c ell interactions to ro bust developmental processes is dis 



cussed as the c ommunity effect 



Gurdon et 



groups of cells (jGreenwald and Rubin 



1993) 



al 



19931 ) and differentiation from equivalent 



When considering the development of a multi-cellular organism, not only a set of cell 
types, but also the number distribution of each of the cell types, has to be suitably deter- 
mined and robust against perturbations during the course of development. The proportion 



(Oviedo et al.. 


2003: 


Rafols et al.. 


2000) 



2000J). In the D. discoideum 



slug, the number ratio of two cell types is kept almost constant. In the hematopoietic sys- 
tem of mammals approximately ten different cell types are generated from a hematopoietic 



stem cell, and their growth and differentiation are regulated to keep the number distri- 
bution of each cell to achieve homeostasis of the hematopoietic system. In this case, in 
addition to the proportion regulation, the absolute size of stem cells is also important 
because all the hematopoietic cells will ultimately die out without their existence. Indeed, 
regulation of the numbers of each cell type is rather common in multicellular organisms. 
As the distribution of each cell type is a property of an ensemble of cells, cell-cell interac- 
tions should be essential for such regulation. 



There are several theoretical studies discussing the importance of cell-cell interactions. 
By considering an ensemble of cells with intra-cellular gene tic (or chemical) ne t works 



and intercellular inte ractions, synchronization of oscillation 



McM 



1995 



lien et al 



2002) or dynamical clusterings 



Kaneko and Yomol . 



19971 : 



Garcia-Oialvo et al 



Kane 



io and Yomo 



Furusawa and Kaneko 



199 



1994; 



2004 



Ullner et al 



Vlizugu chi and Sano , 



2007 



Koseska et al 



20071 ) are observed. Cell states distinguishable from those of a single-cellular dynamics are 
generated, providing a basis for functional differentiation for multicellularity. The preser- 
vation of t he proportion of diffe r ent ce ll types is realized by taking advantage of Turing 



instability (jMizuguchi and Sano 



19951 ). while the robustness i n the number distribution 



of different cell types i s discovered in reaction network models (IKaneko and Yomol . 



Furusawa and Kaneko 



1998 



Kaneko and Yomo 



1994; 



19991 ). Nevertheless, regulatory mecha- 



nisms for cell type populations are not elucidated in terms of dynamical systems because 
of the high dimensionality of the models. 



In the present paper, we propose a regulatory mechanism of cell differentiation based 
on dynamical systems theory by taking simple cell models with biological gene regulation 
dynamics. Specifically, we study how cell states are differentiated with the change in the 
total cell number following cell-cell interactions. By incorporating different interaction 
kinetics, we show how simple functional modules generate specific cellular behaviors such 
as a cell fate switching, size regulation of each cell type, and preservation of the number 
ratio of each cell type. 



3 



The present paper is organized as follows. In section 2, we introduce an interact- 
ing multicellular model which is further analysed in the present paper. Each cell has a 
simple functional module of genes, and its expression dynamics is modulated by the in- 
teractions with other cells. In sections 3 to 5, we consider several different intercellular 
interactions, respectively. Although possible cell states are generated by an intracellular 
functional module, selection of one of these possible states or establishment of specific 
number distributions of cell states is realized depending on the manner of intercellular 
interactions. In section 6, we extend our theoretical scheme to discuss the distribution 



of cell types in cell dif f erent i ation models studied so far (jKaneko and Yomo 



Furusawa and Kaneko 



1998 



1997 



1999 



200ll ). Although these models have a complex intra-cellular 
reaction network, we show that the same logic can be applied to explain the cell differen- 
tiation observed in these models. In section 7, we summarize our results and discuss their 



biological relevance and future directions. 



2 Framework of the Model 

Here we introduce a basic model of interacting cells with intracellular gene expres- 
sion dynamics. Consider N cells with identical genes which interact through a common 
medium. The internal state of i-th cell is represented by the expression pattern of m 
genes, as ui = (uj, . . . , u™) T . The medium under which cells are placed is represented by 
concentrations of n diffuse signals given by v = (v\, . . . ,v n ) T . As the simplest case, we 
discard the spatial configuration of cells so that each cell interacts with all the other cells 
via common signal chemicals v. Each intracellular gene expression dynamics is modulated 
by these signal molecules, which give interactions with other cells. 



For the sake of simplicity, we mostly examine the dynamics of single gene expression, in 
which the state of the i-th cell is expressed by only one variable, Uj, and the intercellular 
interaction is mediated by only one global diffusive signal, v. By using the standard 
kinetics of gene expression, {u^} and v are chosen to obey the following equation, 
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f( UiV ) = -( ?™ U u\ + A \ for i = i n, (1) 

JK ' t \K« + v a (t) + uf{t) w 7 

g(iti,...,ujv, , y)- (2) 

Here gene itj activates its own expression by feedback, while the signal v has an inhibitory 
effect on the expression of the gene ttj. Generally, the signal v is released by each cell 
depending on its gene expression level and the signal abundances at that moment. We 
adopt Hill- type kinetics for self activation of the gene u\. The parameter a denotes the 
Hill coefficient, i.e., the cooperativity of its kinetics, while K u is the threshold for the 
activation of gene Ui, and A u is the activation rate of U{ by other molecules in the cell. 
The parameter r is a time constant of the expression dynamics of m normalized by that 
of the signal v. In the present paper, we consider that the timescale of m is much slower 
than that of v, so that only fixed point solutions are realized. The assumption on the time 
scale is biologically reasonable because the gene expression process requires a much longer 
time than simple catalytic reactions. For numerical simulations, we use the following pa- 
rameter values; K u = 0.1, A u = 0.04, a = 2.0, r = 10.0. Note that the following results 
are qualitatively invariant as long as the Hill-coefficient a is larger than unity. 

Before studying the dynamics of a population of interacting cells, we first survey the 
single intracellular dynamics Eq. (pTJ) with v given as a constant control parameter. As is 
shown straightforwardly, the equation has a fixed point solution which exhibits two saddle- 
node bifurcations with the change in v (Fig. [1]). We denote these bifurcation points as 
v = v\ and v = v\, and call the upper branch of the stable state as itm (or cell state 1) 
that is stable at v < v%, and the other lower branch as um\ (or cell state 2) that is stable 
at v > v*. In the parameter region v* < v < the bistability of tin) and U( 2 ) is sustained. 

As shown in Fig. [H the only possible stationary states of each cell are = urn or 
Ui = U( 2 ). Depending on the value of v and also on the initial condition of Ui, each of the 



duj(t) 

dt 
dv(t) 
dt 
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two solutions are selected. The question we address is as follows: how are these states 
selected and what determines a possible range in the number distribution of the two states 
when intercellular interactions through v are taken into account. In the following sections, 
we analyze three models with different types of the function g(u\, . . . ,un,v) to study how 
the differences in the kinetics of v lead to different types of regulation in the number 
distribution of cell types. 

3 Model I: Cell Fate Determination by Total Cell Number 

As a first example of interacting cells, we adopt a model in which each cell simply 
emits the signal v with the same rate. The kinetics of v obeys the following equation, 



model I 

= gi(«i, ...,u N ,v)=^2ci- v(t) = ciN - v{t), (3) 

i=l 

while the kinetics of {ui} obey Eq. ([1]). We are interested in the behavior of the station- 
ary state as a function of the total cell number N. The stationary state solution of an 
ensemble of cells is generally obtained by the following procedure. First, we regard the 
signal v as a fixed parameter, not a variable, and obtain the solution Ui as a function of 
v, as already described in the previous section. Next, we write down v as a function of 
N and {u^} so that the self-consistent solution of the coupled equation is obtained, from 
which we analyze the dependence of the solution on the total cell number. 

The stationary state is simply obtained by dui/dt = and dv/dt = 0. In the present 
case, the solution v is independent of {itj}, and depends only on N, which leads to 

f(ui,v)=0, v = dN. (4) 

The solution curve f(ui,v(N)) = is shown in Fig. [5J and the numerical result of the 
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ratio of the number of each cell type to the total cell number is shown in Fig. [3J Here 
we define a single-cluster of an ensemble of cells as a state in which all the cells take the 
same stationary states, i.e., 



^ = u {k) (k = 1, or 2) for i = 1, 2, . . . , N, (5) 

and a two-cluster state as that in which two cell types with u = U(\) and u = up) coexist, 
so that 

for % = 1, ...,iV(i), 

(6) 

u (2) for i = N {1) + 1, . . . ,iV (1) + 2V~ (2) (= TV). 
Here iV/^ and -/V( 2 ) denote the number of the cells with u = ury) and u = up) , respectively. 

When the cell number N is lower than a threshold (= v\_/c\), the single-cluster 
state of tin) is realized, while for larger than a threshold iV| (= v^/ci), the single-cluster 
state of U( 2 ) is realized, irrespectively of the initial cell state. Only within the range of 
Ni < N < N£ are two-cluster states of um and up) possible, where any population ratio 
of the cell types with U(x) to up) can be realized depending on the initial condition. Cell 
types switch between um) and up) simply by the total cell number, and the signal v works 
as a population size detector. 



4 Model II: Diversification from Single State, and Size Reg- 
ulation of Specific Cell Type 

Next, we consider the case in which the signal induction depends on the expression 
level of We will show that the cells are differentiated into two types over a wide range of 
the total cell number N, and that the number of type 1 cells remains at a same level herein. 
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The kinetics of the signal v in model II is represented as follows, 



model II 



^ = g 2 ( Ul) . . . , u N , v) = c 2 g - v(t). (7) 

We here adopt Hill- type kinetics for the induction of the signal v by 114 , where (3 is the Hill 
coefficient, representing the cooperativity in the induction, and K v denotes the threshold 
value for the signal induction. The parameter c 2 gives the release rate of v from each cell. 



Dependence of the stationary states on the total cell number is shown in Fig. 01 For 
a small N, all the cells always fall on a single-cluster state of %). As N gets larger, the 
bifurcation to a two-cluster state occurs, where the cells take either or u^)- Here, the 
single-cluster state of ttm (u^)) is realized only at a small (large) number of cells, respec- 
tively, so that there is a gap in the total number of cells between the two single-cluster 
states. The two-cluster state exists within this gap. 



To understand the observed dependence of the clustering behavior on the cell number, 
we first consider the stability of a single-cluster state. From dui/dt = 0, dv/dt = 0, and 
Ui = ut k \ (k = 1, or 2) for i = 1, . . . , N, we get 

/(«(*),*) = 0, v = c 2 N W . (8) 

Kg + u p {k) 

By solving the above equations self-consistently, the solution curve of u is obtained as 
a function of the total cell number iV (Fig. [5]). For N < Nf, a single-cluster state of 
is always stable. When the cell number increases beyond N*, this single-cluster state 
becomes unstable, while for much larger N such that N > iV|, the single-cluster state 
becomes stable again, where the cell state is U( 2 ) (Fig. E|) . The threshold and N% are 
given by JVJ = v* 2 (kS + uLfa))/ fauUv})) and iV| = v\{K^ + «f 2) K))/(oi«f 2) K)), 
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respectively. 



Next, consider the condition for the existence of a two-cluster state. Because the sta- 
bility of a cell state is determined by the amount of v, the condition for the existence of a 
two-cluster state is given by v*[ < v < v 2 ■ Accordingly, considering v as a function of Nn^ 
and N, a two-cluster state is possible if Nn\ satisfies v{ < v(N^, N) < u|. Note that v sat- 
isfies dv{Nn\, N)/dN(i~\ > 0. Thus, the range of the cell number N in which a two-cluster 
state exists is given by N± < N < TVJ, where = v\{K^ + uf^(v*)) /(c2U^(v*)) and 
N 2 = vl(K$ + UpJv^/fau^JvZ)), respectively. These threshold sizes satisfy N* < N* 
and N 2 < N 2 , so that only two-cluster states are stable for N satisfying iV* < N < N%. 



Because the number of each cell type in these two-cluster states has to satisfy the 
above condition, the range of possible numbers of two cell types is limited, depending 
on the total number of cells. The number of cell type 1 (iVm) from a variety of initial 
conditions is plotted as a function of N in Fig. As N is increased beyond N£, TVm 
decreases linearly with N, with a rather small slope, over a wide range of N, up to N 2 . 
Within this range the value of JVi does not change so much. 

To understand this behavior we obtain the dependency of Nry on v and N. In a 
two-cluster state (N(u,N( 2 \), v is expressed by the contribution from the cell types 1 and 
2. Thus, iVn) is written as 

N {1) {N,v) = -A(v)N + B(v), (9) 



u {2) {v)/{K?,+u[ 2) {v)) 
(«)/(*? +«£,(«)) - tf 2) (v)/(K£+u^(v)) 



u 

B{V) = ^{tijj {v)/{Kl + („)) - {v)/(Ki + vf* m («))} ' 
Here, we note that un\ and U( 2 ) are determined self-consistently as functions of v, and 
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that A(v) > and B(v) > 0. For the existence of a two-cluster state, v has to satisfy 
v* < v < u|, that is, N^(N,v*) < N^(N,v) < (N, t>|) for each N. By inserting Eq. 
Q into this expression, it is shown that Nm(N, vf) and N^(N, v 2 ), i.e., the lower and 
upper bounds of ATm, decay linearly with N, with the slope of A(v\) and A(v 2 )- In fact, 
a linear decrease in Nru with the increase in A~ is clearly discernible in Fig. [6l 

Next, we evaluate the value of the slope A{v). Eq. Q is written as A(v) = {(up) I 'Kvf + 
(ti( 2 )/^(i))^}/{l — (u(2)/ u (i))^}- If ""(2) ^ ""(l) an d u (2) ^ ^ are satisfied, that is the 
case for the parameters used in Fig. [6l A(v ) is much smaller than unity. As a result, the 
decrease in Nt±\ with A r is slow, and is sustained at a same level over a wide range 
of A 7 ", satisfying N (1) (vt) < N (1) (v) < N {1] (v* 2 ) (Fig. E}. 

By increasing the Hill-coefficient (3, A(v) becomes much smaller than unity which 
asymptotically go to zero, even if the value of U( 2 ) is the same level as tin) or K v as is 
shown in Fig. [TJ Note that the conditions U( 2 ) < U(\) and U( 2 ) < K v have to be satisfied. 
The value of the slope A(v ) shows an exponential decrease with (5. Hence, is sustained 
at an almost constant level and the population size regulation of cell type 1 is realized 
with a sufficiently large (3. 

5 Model HI: Proportion Preservation of Two Cell Types 

For precise body plan or for tissue homeostasis, proportion regulation of the number 
of each cell type is required. The fraction of each cell type has to be sustained at a certain 
range, against the change in the total number of cells. Here we modify the kinetics of v 
in the previous model II to seek for the possibility of the proportion regulation. With this 
modification, we will show that the population fraction of the two types of cells is kept at 
a certain level against the change of N. 

Here, the kinetics of v is modified as follows, 
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model m 



dv(t) 

—fij- = g3(«l, ...,Un,V) 

N v?(t) N k 13 

^w^-^ v{t) EwTW)- v(t) (12) 

The modification to model II is just an addition of the second term in Eq. (|12[) . In other 
words, each cell in this model also contributes to the degradation of the signal v. 



As in the previous model, the cellular states fall on stationary states, and the bifur- 
cation of the stationary state from a single-cluster to two-cluster states are observed with 
the increase in A?" (Fig. [8j). Here, we first note that the two-cluster state remains stable 
over a wide range of N. Indeed, non-zero exists so that u* < v(N^,N) < is 
satisfied even for sufficiently large N. 



Next, we study the population distribution of two cell types. As shown in Fig. [9l the 
ratio N^/N stays at a constant level against the change of N. In the same way as in the 
previous section, the dependency of Nn\ on v and N for a two-cluster state is written as, 

^ = a {v) + M, (13 ) 



A{y) 
B{v) 



(14) 



=„i {«f 1} W/iKg+uf^ («))-«f 2) (v)/(Kg+uP 2) (v))}+c v2 Kgv{ l/(Kg+uf 2) (v))-l/{Ki +uf 1} («))} ' 



(15) 



Here, B(v) > is always satisfied. Because v satisfies v\ < v < v\ for the existence of a 
two-cluster state, N^/N is within the range (A(v*)+b(vD/n) < N wj^' v * > < [A(v^)+b(v* 2 )/n) 
for each N. As a result, when A" is sufficiently large, the possible range of N^/N is given 
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by 

m) < ^ < (i6) 

Prom the above expression of A(v), if the condition (u^/^fi) (^2)) < Cvi/(cv2Kv) < 
(u|/u/L ( u i)) * s satisfied, A(u) is within < A(v) < 1. This is the case for the parameter 
values in Fig. [9j Thus, the cell type ratio of a two-cluster state has to be within the range 
given by Eq. (|16p . so that its ratio is insensitive to the change of the total number of 
cells. In addition, by increasing the Hill-coefficient P, the range given by Eq. (|16p gets 
narrower. Thus, the ratio Nru/N is more accurately regulated. As (3 goes to infinity the 
range approaches its minimum, where the boundary is given by A^v) = v/{c v \/c V 2 + v). 

Note that A(v) here is positive and is not necessarily small, in contrast to A{v) in Eq. 
(J9]) for the model II. Inclusion of the second term in Eq. (|12|) allows for this behavior, and 
the proportion regulation of cell types is achieved over a wide range of cells. 



6 Cell Differentiation Model with Random Network 



Here, we briefly discuss a general situation of cell differentiation models with in- 
tracellular dynamics and intercellular interactions with more genes (chemical species). 
As an e xample, We use th e cell differentiation models of Kaneko-Yomo o r Furusawa- 



Kaneko (jKaneko and Yomol . 



1997 



1999 



Furusawa and Kaneko 



1998 



2001 



Here, we 



aim at demonstrating that the regulative behavior of cell differentiation in the previous 
sections generally works, which at the same time may provide a possible explanation for 
differentiation phenomena observed in their m odels. For the following ana lysis, we use 



one of the models (FK model) introduced in ([Furusawa and Kanekol . 
straightforwardly extended to other models. 



2001 



, while it is 



In the FK model each cell has intracellular metabolic dynamics, and grows by uptake 
of the nutrients in the medium, and divides when the abundances of chemicals in the 
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cell goes beyond some threshold. Accordingly, the total cell number N is also a time- 
dependent variable. As the cells share the same medium, they interact with other cells 
through uptake from the medium and exchange of chemicals with it. 



The state of cell I is expressed by P different metabolites, = (x\ , . . . ,Xp) T , and 
the nutrients are X = (X±, . . . , Xq) t , (Q < P). The dynamics of the i-th metabolite in 



cell I is given as follows 



(17) 



dt 



A change in the concentration of the i-th nutrient in the medium with the volume V is 
given by, 



Si is the external source of the nutrient, D env is the diffusion constant between the nutri- 
ent reservoir and the medium, and D is that across the cell membrane. Each cell grows 
through uptake of nutrients and changing them to other metabolites by Eq. (|17|) . As 
the cells share the same medium, they interact with each other through competition for 
nutrients. Here we confine our consideration only to the behavior of nutrients {X{\ in 
the stationary states for fixed N, and to obtain the behavior of the stationary states as a 
function of N. 

Because the stationary state satisfies the condition dxf' /dt = and dXi/dt = 0, 



dXi(t) 



D cnv (S l -X i (t))--Y,(X i (t) 




(18) 



771=1 




(19) 



(20) 



»i=i 



From Eq. (|19|) . possible stationary states of each cell, i 



.e., stationary solutions of xft' 



are 



obtained as a function of X. Next, we describe how X varies with N. As in the previous 
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sections, we assume that the cell population takes an M-cluster state in the stationary 
state for a given N. By solving Eq. (|20p for JQ, one obtains 



X: 



jV 



1 + 



- 



(21) 



- (k) 

where x\ ' is the cell type k in an M-cluster state, and = N^/N, with as the number 
of type k cells in the population. Xj is represented as a function of N and The stabil- 

ity condition of the M-cluster state of concern is expressed by N and {Rk}, from Eq. (|19j) 
and Eq. (|21|) . Thus, the realization of an M-cluster state depends on the number of cells 



or the ratio of cell types. Regulation of each ce ll type, as observed in (jKaneko and Yomol . 



1997 



1999: 



Furusawa and Kanekol . 



1998 



2001 



is expected accordingly. 



7 Summary and Discussion 

Through the analysis of several models, we see, i) a switch of cell types via an increase 
of the total cell number, and ii) diversification to two cell types. In addition, when the cells 
differentiate to two types, size preservation of a specific cell type or proportion preservation 
of two cell types appears, depending on the interaction form with other cells. These be- 
haviors are explained as a bifurcation of cell states via the intercellular interactions. First, 
possible cell types and u^) are generated by a single positive feedback loop, which 
works as a module for bistability. Secondly, intercellular signal v works as a bifurcation 
parameter, whose abundances determine the actual cell types. This bifurcation parameter 
is a function of the number of each cell type, depending on the intercellular interactions. 
Then, the resulting bifurcation parameter has to be determined self-consistently. This 
constraint restricts the number distribution of the cell types, which gives the mechanism 
of the regulation of the cell differentiation. 



In model I, because the total cell number simply corresponds to the bifurcation param- 
eter of cell states, the switch of the cell types by the total cell number is straightforward. 
In model II and HI, since intercellular couplings change the bifurcation parameter, the 
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transition from the single-cluster state of urn to a two-cluster state occurs by the increase 
in the total cell number. In model II, the cell-type 2 contributes only weakly to the increase 
of v, compared with the cell-type 1. Thus, the amount of v mainly depends on the number 
of the cell-type 1. In contrast, in model HI, the cell-type 2 degrades v. As a result, the 
amount of v depends on the number ratio of two cell-types. 



If a gene expression network shows bistability with a bifurcation structure as in Fig. 
1, cell differentiation is a general consequence when cell-cell couplings are introduced. An 
important point here is that the same intracellular module can be used in several different 
biological contexts by modifying only the intercellular interaction. This is quite useful in 
an evolutionary perspective because new biological functions can be added by incorporat- 
ing new interactions while preserving the intracellular core module. 



Here, we discuss several biological examples that may correspond to our models. First, 
we refer to the cell cycle machinery in Xenopus, where a Cdc2 positive feedback loop makes 
a bifurcation with regards to the amount of cyclin B, and introduces bistability as in Fig. 
1. In this system, an increase in the DNA amount in the early embryo induces the transi- 
tion from the low Tyrl5 phosphorylation state of Cdc2 to the high Tyr l5 phosphorylation 



state , which seems t o cause the mic 



1993 



Hartley et al. 



1996 



Sha et al. 



jlastu la transition in Xenopus (jNovak and Tyson 



20031 ). The induced differentiation may correspond 



to that observed in model I. 



Secondly, as for model II, consider the maintenance of the h e mato poietic stem cells 



2009 ). The stem cells 



in mammals, where osteoblasts work as a stem cell niche (jCalvil . 
compete for some chemical factor representing this niche, and the cells which cannot take 
the factor differentiate to specific hematopoietic lineages. It has been discussed that the 
regulation of the stem cell population size is realized through the competiti on for the fac 



tor, to which responsibility decr eases through the differentiation process (jRadtke et al. 



2004 



Adams and Scadden 



20061 ). Indeed, it is observed that the expression of Notchl, 



15 



which is a candidate for the involveme nt in the nich e -stem interaction, disappears after 
commitment to the lymphoid lineage (jRadtke et all 120041 ). The differentiation in the 
hematopoietic system may correspond to that studied in model II. 



Thirdly, an example for model III is given by the proportion regulation of prestalk- 
cell types and prespore-cell types in the Dictyostelium slug. Differentiation to prespore 
cells is induced by cAMP, and the cell state is mai ntained by a positive- feedbac k loop 



of prespore cell spec i fic ad enylyl cyclase G activity (jHopper et al 



Alvarez-Curto et al 



1993; 



Williams 



2006; 



20071 ) . On the other hand, differentiation-inducing factor-1 (DIF- 



1) is necessary for the differentiation from a prespore-cell to a prest alk-cell (a t least 



for the differentiation to p stO which is a subtype of the prestalk-cell) ([Williams 



Kay and Thompson 



2001 



2006 



). As an intercellular interaction, this DIF-1 is produced by 
prespore-cells, and are degraded by prestalk-cells. This cell-type specific induction/destruction 
of DIF-1 is responsible for the proportion preservation as studied in model HI. 



Although we confine our analysis to a system with only fixed point solutions, oscilla- 
tory and other dynamical behaviors are often observed in biological systems. The analysis 
we introduced here is also applicable to such cases, as long as there are bifurcations of 
attractors with the change in relevant chemical concentrations which are influenced by 
cell-cell interactions. On the other hand, oscillatory behaviors may bring about richer bi- 
furcations, as well as clustering of cells with regards to the oscillation phase or amplitude, 



as has been discussed in mode 



ter actions ( Kaneko and Yomo 



s wit 



i infra-cellular oscil 



1994, 



1997 



Koseska et al 



atory 



2007 



dynamics and cell-c ell in- 



Ullner et al 



2007J). The 



study of possible forms on differentiations and regulations in such dynamical systems will 
be important in future. In multicellular systems, cells behave in coordination by taking 
advantage of communication with other cells. Such collective behavior is a result of inter- 
acting systems with intra-cellular gene expression dynamics. The present self-consistent 
determination of bifurcation parameters through cell-cell interactions will be essential to 
understand organization in multicellularity. 
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Captions 

Figure 1: 

The value u of the fixed point solution as a function of the signal concentration v in 
our model. Solid line indicates the stable solution, while the dotted line indicates the 
unstable one. 

Figure 2: 

The stationary states of Ui in model I are plotted against the total cell number N. At 
the interval N* < N < iV|, two different cell states coexist. The parameter value c\ is set 
at 0.005. 

Figure 3: 

The ratio of the number of each cell type (x for TVm and □ for Nm) plotted against 
the total cell number N, for model I. The initial values of m are chosen randomly from 
the interval of ui € [0, 1]. The parameter value is c\ = 0.005. 

Figure 4: 

The fixed point values of Ui in model II are plotted against the total cell number N . At 
each N, 100 initial conditions are chosen. The expression levels of Ui for a single cluster 
(+) and two-cluster solutions (o) are plotted as a function of N . The value for two-cluster 
solutions is the average over initial conditions. The parameter values are set at K v = 2.0, 
(3 = 2.0, c 2 = 0.1. 

Figure 5: 

The stationary state of a single-cluster solution for model II. Solid line indicates Ui 
of the stable fixed solution, while the broken line denotes that of the unstable one. The 
parameters are K v = 2.0, (3 = 2.0, c 2 = 0.1. 

Figure 6: 
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The number of cell type 1 (x) is plotted against the total cell number in model II. 
The initial condition of itj is chosen randomly from the interval Ui £ [0, 1]. Solid and bro- 
ken lines indicate N^(N,vl) = -A(vl)N + B(x>l) and iV (1) (iV,u*) = -A{v%)N + B(v%), 
respectively, where = 0.011, B{v\) = 24, A{v* 2 ) = 0.0099, and B(v* 2 ) = 93. The 

parameters are K v = 2.0, (3 = 2.0, C2 = 0.1. 

Figure 7: 

The slope A(f3\v) is plotted as a function of f3. Here, A((3;v) for two different values 
of v, i.e., v* and v\ are plotted, which agree within the resolution of the plot in the figure. 
The parameter values are K v = 2.0, C2 = 0.1. 

Figure 8: 

The fixed point solutions of model IH plotted against the total cell number N. At each 
N, 100 initial conditions are chosen. The expression levels of Ui for a single cluster (+) 
and two-cluster solutions (o) are plotted as a function of N. The value for two-cluster so- 
lutions is the average over initial conditions. The parameter values are K v = 0.2, (3 = 2.0, 
Cvi = c V 2 = 0.005. 

Figure 9: 

The ratio of the number cell type 1 to the total cell number N is plotted against N 
for model HI. The initial condition of ui is chosen randomly from the interval of Uj € [0, 1]. 
Solid and broken lines indicate N {1) (N,vl)/N = A{v{) + B{v\)/N and N {l) {N ,v* 2 )/N = 
A(v%) + B(v* 2 )/N, respectively, where A(v*) = 0.16, B(vf) = 69, A(v* 2 ) = 0.36, and 
B(v 2 ) = 86. The parameter values are K v = 0.2, (3 = 2.0, c v \ = c V 2 = 0.005. 
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